!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!
!
C
C  /* Deck cc_soldrv */
C--------------------------------------------------------------------
      SUBROUTINE CC_SOLDRV(LIST,NSTAT,ORDER,ISIDE,APROXR12,
     &                     ISYMS, ISTAT, EIGVAL, LPROJ,
     &                     ISYMO, LABEL, FREQS, ICAU,
     &                     ISYOF, NVEC,  MAXVEC, 
     &                     WORK, LWORK)
C--------------------------------------------------------------------
C
C     Purpose: Direct solution of Coupled Cluster response equations.
C
C              LIST  -- list name, used to get the gradient vector and
C                       to read/write result to file
C
C              NSTAT -- nb. of involved excited states
C              ORDER -- response order (n.b. of Labels,frequencies...)
C
C              ISIDE = 1 : solve right response equations (A-w)R = X
C              ISIDE =-1 : solve left response equations  L(A+w) = X
C
C              ISYMS   -- symmetry of involved excitated states
C              ISTAT   -- indices of excited states
C              EIGVAL  -- excitation energies
C
C              ISYMO   -- operator symmetries
C              LABEL   -- operator labels
C              FREQS   -- associated frequencies
C              ICAU    -- cauchy order
C
C              ISYOFF  -- symmetry offsets, as generated by LSTSORT
C              NVEC    -- number of response equations to solve
C              MAXVEC  -- first dimension in allocation of ISYMS
C                         ISTAT,EIGVAL,ISYMO,LABEL,FREQS,ICAU
C
C              WORK/LWORK -- work space, length of work space
C            
C
C  Christof Haettig, Maj 1997, based on Ove's version of CC_SOLRSP
C  Christof Haettig, Jul 1997, excited state resp. eq. implemented
C  Christof Haettig, Oct 1997, batching for Cauchy eq. implemented
C  Christof Haettig, Feb 1998, second-order Cauchy eq. implemented
C  Sonia Coriani,    Mar 2000, added PL1 to ELn/ERn special batching
C  Christof Haettig, Apr 2002, small modifications for CC3
C----------------------------------------------------------------------
C
#include "implicit.h"
#include "priunit.h"
#include "dummy.h"
#include "maxorb.h"
#include "iratdef.h"
#include "cclr.h"
#include "ccorb.h"
#include "ccsdsym.h"
#include "ccfro.h"
#include "ccsdio.h"
#include "ccsdinp.h"
#include "leinf.h"
Cholesky
#include "ccdeco.h"
Cholesky

C
      LOGICAL LOCDBG
      PARAMETER (LOCDBG = .FALSE.) 
C
      LOGICAL TRIPLET
      PARAMETER (TRIPLET= .FALSE.) 
C
      CHARACTER LIST*(*), APROXR12*(*)

      INTEGER NSTAT, ORDER, MAXVEC, NVEC, ISIDE, NSIMLEBAK
      INTEGER ISYOF(8)

      LOGICAL LPROJ(MAXVEC)
      CHARACTER*8 LABEL(MAXVEC,ORDER)
      INTEGER ISYMS(MAXVEC,NSTAT), ISTAT(MAXVEC,NSTAT)
      INTEGER ISYMO(MAXVEC,ORDER), ICAU(MAXVEC,ORDER)

      DOUBLE PRECISION FREQS(MAXVEC,ORDER), EIGVAL(MAXVEC,NSTAT)
      DOUBLE PRECISION SIGN
      DOUBLE PRECISION ZERO, ONE
      PARAMETER (ZERO = 0.0d0, ONE = 1.0d0)

      CHARACTER*8 FC1AM,FC2AM,FRHO1,FRHO2,FRHO12,FC12AM,FS12AM,FS2AM
      PARAMETER (FC1AM ='CCR_C1AM',FC2AM ='CCR_C2AM')
      PARAMETER (FRHO1 ='CCR_RHO1',FRHO2 ='CCR_RHO2')
      PARAMETER (FRHO12='CCR_RH12',FC12AM='CCR_C12M',FS12AM='CCR_S12M')
      PARAMETER (FS2AM='CCR_S2DM')

      LOGICAL LINQCC
      PARAMETER (LINQCC = .TRUE.) 

      DIMENSION WORK(LWORK)
      CHARACTER MODEL*10
      LOGICAL LTAKE, LPROJECT, LREDS

Cholesky
      PARAMETER (NCHCHK = 2)
      CHARACTER*3 CHOCHK(NCHCHK)
      DATA CHOCHK /'L0 ','R1 '/
Cholesky

C
      CALL QENTER('CC_SOLDRV')
      LUFC1  = -1
      LUFC2  = -1
      LUFC12 = -1
      LUFR1  = -1
      LUFR2  = -1
      LUFR12 = -1
      LUFS12 = -1
      LUFS2  = -1
C
      LREDS = CCR12
C
      IF (LOCDBG) THEN
        WRITE (LUPRI,*) 'LIST  :',LIST
        WRITE (LUPRI,*) 'NSTAT :',NSTAT
        WRITE (LUPRI,*) 'ORDER :',ORDER
        WRITE (LUPRI,*) 'ISIDE :',ISIDE
        WRITE (LUPRI,*) 'NVEC  :',NVEC
        WRITE (LUPRI,*) 'MAXVEC:',MAXVEC
        WRITE (LUPRI,*) 'CCR12,LREDS:',CCR12,LREDS
        WRITE (LUPRI,*) 'ISYMS :',((ISYMS(I,J),J=1,NSTAT),I=1,NVEC)
        WRITE (LUPRI,*) 'ISTAT :',((ISTAT(I,J),J=1,NSTAT),I=1,NVEC)
        WRITE (LUPRI,*) 'EIGVAL:',((EIGVAL(I,J),J=1,NSTAT),I=1,NVEC)
        WRITE (LUPRI,*) 'ISYMO :',((ISYMO(I,J),J=1,ORDER),I=1,NVEC)
        WRITE (LUPRI,*) 'LABEL :',((LABEL(I,J),J=1,ORDER),I=1,NVEC)
        IF (LIST(1:1).NE.'C' .AND. LIST(2:2).NE.'C' .AND.
     &      LIST(1:2).NE.'M1'.AND. LIST(1:2).NE.'N2' ) THEN
          WRITE (LUPRI,*) 'FREQS :',((FREQS(I,J),J=1,ORDER),I=1,NVEC)
        END IF
        IF (LIST(1:1).EQ.'C' .OR.  LIST(2:2).EQ.'C') THEN
          WRITE (LUPRI,*) 'ICAU  :',((ICAU(I,J),J=1,ORDER),I=1,NVEC)
        END IF
        IF (LIST(1:2).EQ.'ER' .OR. LIST(1:2).EQ.'EL') THEN
          WRITE (LUPRI,*) 'LPROJ :',(LPROJ(I),I=1,NVEC)
        END IF
      END IF

      IF (NVEC.EQ.0) THEN
         CALL QEXIT('CC_SOLDRV')
         RETURN
      END IF

C
C-------------------------------------
C     Initialize Leinfi and variables.
C-------------------------------------
C
      IPRLE = IPRINT 
      MODEL = 'xxxxxxxxxx'
      IF (CCR12) THEN
        IF (CCS)             MODEL = 'CCS       '
        IF (CC2 .AND. CCR12) MODEL = 'CC2-R12   '
        IF (CCSD.AND. CCR12) MODEL = 'CCSD-R12  '
        IF (CC3 .AND. CCR12) MODEL = 'CC3-R12   '
      ELSE
        IF (CCS)  MODEL = 'CCS       '
        IF (CC2)  MODEL = 'CC2       '
        IF (CCD)  MODEL = 'CCD       '
        IF (RCCD) MODEL = 'RCCD=RPA  '
        IF (DRCCD) THEN
           IF (SOSEX) THEN
                  MODEL = 'SOSEX     '
           ELSE
                  MODEL = 'DRCCD     '
           END IF
        END IF
        IF (CCSD) MODEL = 'CCSD      '
        IF (CC3)  MODEL = 'CC3       '
      END IF
C
C      IF (.NOT.(CCS.OR.CC2.OR.CCSD.OR.CC3)) THEN
      IF (.NOT.(CCS.OR.CC2.OR.CCD.OR.CCSD.OR.RCCD.OR.DRCCD.OR.CC3)) THEN
         WRITE(LUPRI,'(A)') ' CC_SOLDRV: Do not want to calculate '
     *           //'anything else than CCS, CC2 '
     *           //'CCD, RCCD, DRCCD, CCSD and CC3 properties '
         CALL QUIT('Model not CCS,CC2,(RD)CCD,CCSD,CC3 in CC_SOLRSP')
      ENDIF
 
C     Cholesky Check Section.
C     =======================
C-tbp: For Cholesky, most of the equations are not yet implemented,
C      and CC2 is the only possible model! Triplet? nope...

      IF (CHOINT) THEN
         IF (.NOT. CC2) THEN
            WRITE(LUPRI,'(A,A,A)')
     &      ' CC_SOLDRV: Cholesky decomposed integrals not implemented',
     &      ' for ',MODEL
            CALL QUIT('Cholesky '//MODEL
     &                //' not implemented in CC_SOLDRV')
         ENDIF
         IF (TRIPLET) THEN
            WRITE(LUPRI,'(A,A)')
     &      ' CC_SOLDRV: Cholesky decomposed integrals not implemented',
     &      ' for triplets'
            CALL QUIT('Triplet Cholesky not implemented in CC_SOLDRV')
         ENDIF
         LFOUND = 0
         DO ICHECK = 1,NCHCHK
            IF (LIST(1:3) .EQ. CHOCHK(ICHECK)) THEN
               LFOUND = 1
               GO TO 111
            ENDIF
         ENDDO
  111    CONTINUE
         IF (LFOUND .EQ. 0) THEN
            WRITE(LUPRI,'(A,A,A)')
     &      ' CC_SOLDRV: Cholesky decomposed integrals not implemented',
     &      ' for ',LIST
            CALL QUIT('Cholesky '//LIST(1:3)
     &                //' eqs. not implemented in CC_SOLDRV')
         ENDIF
      ENDIF

C---------------------------------------------------------------------
C        For CC3: switch batching OFF 
C        (presently we can not solve CC3 equations simultaneously,
C         because then we are not allowed to share the same reduced
C         space for equations with different frequencies...)
C---------------------------------------------------------------------
      NSIMLEBAK = NSIMLE
      IF (CCSDT) NSIMLE = 1

C-------------------------------------------------------------
C     Start loop structure - first loop over symmetry classes.
C-------------------------------------------------------------
 
      DO ISYM = 1, NSYM
 
C-----------------------------
C        Initialize Variables.
C-----------------------------
 
         ISYMTR = ISYM
Chol     IF (CCS) THEN
         IF (CCS .OR. (CHOINT .AND. CC2)) THEN
           NCCVAR = NT1AM(ISYMTR)
         ELSE
           NCCVAR = NT1AM(ISYMTR) + NT2AM(ISYMTR)
           IF (CCR12) NCCVAR = NCCVAR + NTR12AM(ISYMTR)
         END IF
         CALL CCLR_LEINFI(TRIPLET)
         THRLE = THRLEQ
      IF (LETYPA .GT.0) THEN
 
C-------------------------------------------------------------------
C        Find out how many equations to solve in this symmetry class
C-------------------------------------------------------------------
         IBOFF = ISYOF(ISYM)
         IF (ISYM .EQ. NSYM) THEN
            IBEND = NVEC
         ELSE
            IBEND = ISYOF(ISYM+1)
         ENDIF
         NEQSYM  = IBEND - IBOFF
 
C---------------------------------------------------------------------
C        Make batching over nr. of simultaneous vectors to be solved:
C        all or take the number from input.
C        For Cholesky CC2 the Jacobian is frequency-dependent. Thus,
C        we solve the equations each frequency. (I.e., the batch
C        length must be the number of unique operators in each sym.).
C---------------------------------------------------------------------
         IF (NEQSYM .EQ. 0 ) THEN
            NSIM = 0
            NBAT = 0
         ELSE
            IF (CHOINT) THEN
               CALL CC_CHONSIM(LIST,ISYM,NEQSYM,NSIMLE,NSIM)
            ELSE
               IF ( NSIMLE .EQ. 0 ) THEN
                  NSIM  = NEQSYM
               ELSE
                  NSIM  = NSIMLE
               END IF
            ENDIF
            NBAT  = (NEQSYM-1)/NSIM + 1
         ENDIF

*---------------------------------------------------------------------
*        special batching for excited state response vectors:
*        and projected Lagragian multipliers
*---------------------------------------------------------------------
         IF (LIST(1:2).EQ.'ER' .OR. LIST(1:2).EQ.'EL'
     &                         .OR. LIST(1:3).EQ.'PL1') THEN
           IREND = IBOFF
           NBAT  = 0
           DO WHILE (IREND.LT.IBOFF+NEQSYM)
               NBAT  = NBAT + 1
               IRST  = IREND+1
               NSIMR = 1
               IREND = IRST
               IF ( IREND .LT. (IBOFF+NEQSYM) ) THEN
                 LTAKE = ( ( ISTAT(IREND+1,1).EQ.ISTAT(IRST,1)
     &                        .AND. (LPROJ(IREND+1).EQV.LPROJ(IRST)) )
     &                    .OR. ((LPROJ(IREND+1).EQV..FALSE.)
     &                           .AND. (LPROJ(IRST).EQV..FALSE.)) ) 
               END IF
               DO WHILE (  IREND.LT.(IBOFF+NEQSYM) .AND.
     &                     NSIMR.LT.NSIM  .AND. LTAKE ) 
                 IREND = IREND + 1
                 NSIMR = NSIMR + 1
                 IF ( IREND .LT. (IBOFF+NEQSYM) ) THEN
                   LTAKE = ( ( ISTAT(IREND+1,1).EQ.ISTAT(IRST,1)
     &                         .AND. (LPROJ(IREND+1).EQV.LPROJ(IRST)) )
     &                      .OR. ((LPROJ(IREND+1).EQV..FALSE.)
     &                             .AND. (LPROJ(IRST).EQV..FALSE.)) ) 
                 END IF
               END DO
           END DO
         END IF

*---------------------------------------------------------------------
*        special batching for Cauchy vectors:
*---------------------------------------------------------------------
         IF ( (LIST(1:2).EQ.'RC'.AND. .NOT.NEWCAU) 
     &        .OR. LIST(1:2).EQ.'LC'
     &        .OR. LIST(1:2).EQ.'CR'
     &        .OR. LIST(1:2).EQ.'CL'               ) THEN
           IREND = IBOFF
           NBAT  = 0
           DO WHILE (IREND.LT.IBOFF+NEQSYM)
               NBAT  = NBAT + 1
               IRST  = IREND+1
               NSIMR = 1
               IREND = IRST
               IF ( IREND .LT. (IBOFF+NEQSYM) ) THEN
C                LTAKE = .TRUE.
C                DO IORD = 1, ORDER
C                  LTAKE = LTAKE .AND. 
C    &                       (ICAU(IREND+1,IORD) .EQ. ICAU(IRST,IORD))
C                END DO
                 ICAUS = 0
                 ICAUE = 0
                 DO IORD = 1, ORDER
                   ICAUS = ICAUS + ICAU(IRST,IORD)
                   ICAUE = ICAUE + ICAU(IREND+1,IORD)
                 END DO
                 LTAKE = ICAUE .EQ. ICAUS 
               END IF
               DO WHILE (  IREND.LT.(IBOFF+NEQSYM) .AND.
     &                     NSIMR.LT.NSIM  .AND. LTAKE ) 
                 IREND = IREND + 1
                 NSIMR = NSIMR + 1
                 IF ( IREND .LT. (IBOFF+NEQSYM) ) THEN
C                  LTAKE = .TRUE.
C                  DO IORD = 1, ORDER
C                    LTAKE = LTAKE .AND.
C    &                       (ICAU(IREND+1,IORD) .EQ. ICAU(IRST,IORD))
C                  END DO
                   ICAUS = 0
                   ICAUE = 0
                   DO IORD = 1, ORDER
                     ICAUS = ICAUS + ICAU(IRST,IORD)
                     ICAUE = ICAUE + ICAU(IREND+1,IORD)
                   END DO
                   LTAKE = ICAUE .EQ. ICAUS 
                 END IF
               END DO
           END DO
         END IF
*---------------------------------------------------------------------


C
         IF (LOCDBG ) THEN
            WRITE(LUPRI,*) '         '
            WRITE(LUPRI,*) 'CC_SOLRSP-1: NEQSYM = ',NEQSYM
            WRITE(LUPRI,*) 'CC_SOLRSP-1: IBOFF  = ',IBOFF
            WRITE(LUPRI,*) 'CC_SOLRSP-1: IBEND  = ',IBEND
            WRITE(LUPRI,*) 'CC_SOLRSP-1: NSIM   = ',NSIM
            WRITE(LUPRI,*) 'CC_SOLRSP-1: NBAT   = ',NBAT
         ENDIF
C
        IREND = IBOFF
        DO IRB = 1, NBAT
 
C----------------------------------------------------
C           Find start number for this batch on list.
C           and the number of equations.
C----------------------------------------------------

*..................................................................... 
*              special batching for excited state response vectors
*              and projected 1st-order Langrangian multipliers:
*..................................................................... 
            IF ( LIST(1:2).EQ.'ER' .OR. LIST(1:2).EQ.'EL'
     &                             .OR. LIST(1:3).EQ.'PL1' ) THEN
               IRST  = IREND+1
               NSIMR = 1
               IREND = IRST
               IF ( IREND .LT. (IBOFF+NEQSYM) ) THEN
                 LTAKE = ( ( ISTAT(IREND+1,1).EQ.ISTAT(IRST,1)
     &                        .AND. (LPROJ(IREND+1).EQV.LPROJ(IRST)) )
     &                    .OR. ((LPROJ(IREND+1).EQV..FALSE.)
     &                           .AND. (LPROJ(IRST).EQV..FALSE.)) ) 
               END IF
               DO WHILE (  IREND.LT.(IBOFF+NEQSYM) .AND.
     &                     NSIMR.LT.NSIM  .AND. LTAKE ) 
                 IREND = IREND + 1
                 NSIMR = NSIMR + 1
                 IF ( IREND .LT. (IBOFF+NEQSYM) ) THEN
                   LTAKE = ( ( ISTAT(IREND+1,1).EQ.ISTAT(IRST,1)
     &                          .AND. (LPROJ(IREND+1).EQV.LPROJ(IRST)))
     &                      .OR. ((LPROJ(IREND+1).EQV..FALSE.)
     &                             .AND. (LPROJ(IRST).EQV..FALSE.)) ) 
                 END IF
               END DO

*..................................................................... 
*              special batching for Cauchy equations:
*..................................................................... 
            ELSE IF ( (LIST(1:2).EQ.'RC' .AND. .NOT.NEWCAU)
     &                .OR. LIST(1:2).EQ.'LC'
     &                .OR. LIST(1:2).EQ.'CR'
     &                .OR. LIST(1:2).EQ.'CL'                 ) THEN
               IRST  = IREND+1
               NSIMR = 1
               IREND = IRST
               IF ( IREND .LT. (IBOFF+NEQSYM) ) THEN
C                LTAKE = .TRUE.
C                DO IORD = 1, ORDER
C                  LTAKE = LTAKE .AND.
C    &                     (ICAU(IREND+1,IORD) .EQ. ICAU(IRST,IORD))
C                END DO
                 ICAUS = 0
                 ICAUE = 0
                 DO IORD = 1, ORDER
                   ICAUS = ICAUS + ICAU(IRST,IORD)
                   ICAUE = ICAUE + ICAU(IREND+1,IORD)
                 END DO
                 LTAKE = ICAUE .EQ. ICAUS 
               END IF
               DO WHILE (  IREND.LT.(IBOFF+NEQSYM) .AND.
     &                     NSIMR.LT.NSIM  .AND. LTAKE ) 
                 IREND = IREND + 1
                 NSIMR = NSIMR + 1
                 IF ( IREND .LT. (IBOFF+NEQSYM) ) THEN
C                  LTAKE = .TRUE.
C                  DO IORD = 1, ORDER
C                    LTAKE = LTAKE .AND.
C    &                       (ICAU(IREND+1,IORD) .EQ. ICAU(IRST,IORD))
C                  END DO
                   ICAUS = 0
                   ICAUE = 0
                   DO IORD = 1, ORDER
                     ICAUS = ICAUS + ICAU(IRST,IORD)
                     ICAUE = ICAUE + ICAU(IREND+1,IORD)
                   END DO
                   LTAKE = ICAUE .EQ. ICAUS 
                 END IF
               END DO
    
               IF (CC3) THEN
                 IF (LOCDBG) THEN
                   WRITE(LUPRI,*) 'GET CAUCHY RHS VECTORS:' 
                   WRITE(LUPRI,*) 'VECTOR LIST :',LIST(1:3)
                   WRITE(LUPRI,*) 'CAUCHY ORDER:',ICAUS
                 END IF

                 IF (LIST(1:2).EQ.'RC') THEN
                    CALL CCRHSVEC('CO1',LABEL,
     &                            IDUMMY,IDUMMY,DUMMY,
     &                            ISYMO,DUMMY,LDUMMY,
     &                            ICAU,NSIMR,MAXVEC,IRST-1,
     &                            WORK,LWORK)
                 ELSE
                    CALL QUIT('Wrong Cauchy vector list in CC_SOLDRV.')
                 END IF

               END IF

               IF (LOCDBG) THEN
                 WRITE(LUPRI,*) 'SOLVE FOR CAUCHY VECTORS:' 
                 WRITE(LUPRI,*) 'VECTOR LIST :',LIST(1:3)
                 WRITE(LUPRI,*) 'CAUCHY ORDER:',ICAUS
               END IF

*..................................................................... 
*              standard batching for ground state response equations:
*..................................................................... 
            ELSE
               IRST  = IBOFF + (IRB-1)*NSIM + 1
               NSIMR = MIN(NSIM,NEQSYM - (IRB-1)*NSIM)
               IREND = IRST  + NSIMR - 1
            END IF
C
C
            IF (LOCDBG ) THEN
               WRITE(LUPRI,*) 'CC_SOLRSP-1: IRST   = ',IRST
               WRITE(LUPRI,*) 'CC_SOLRSP-1: IREND  = ',IREND
               WRITE(LUPRI,*) 'CC_SOLRSP-1: NSIMR  = ',NSIMR
            ENDIF
C
            IF (IPRINT.GT.2 .OR. LOCDBG) THEN
               WRITE (LUPRI,'(/,1X,A,I3,A)')
     &            'Solving ',NSIMR,' response equations:'
               DO IR = IRST, IREND
                 IF (LIST(1:1).EQ.'C' .OR. LIST(2:2).EQ.'C') THEN
                   WRITE (LUPRI,'(1X,2A,10(A8,I3,I3))') LIST,':',
     &               (LABEL(IR,I),ICAU(IR,I),ISYMO(IR,I),I=1,ORDER)
                 ELSE
                   WRITE (LUPRI,'(1X,2A,10(A8,F10.6,I3))') LIST,':',
     &               (LABEL(IR,I),FREQS(IR,I),ISYMO(IR,I),I=1,ORDER)
                 END IF
               END DO
            ENDIF
            IF (IPRINT.GT.10 .OR. LOCDBG)
     *            WRITE(LUPRI,*) 'CC_SOLRSP-1: Workspace:',LWORK
C
            CALL FLSHFO(LUPRI)
C
C------------------------------------------------------------------
C           Allocation of work space.
C           NB: The gradient vectors is not in memory at this time.
C------------------------------------------------------------------
C
            KIPLAC = 1
            KREDH  = KIPLAC + MAXRED
            KREDGD = KREDH  + MAXRED*MAXRED
            KEIVAL = KREDGD + MAXRED*NSIMR
            KSOLEQ = KEIVAL + MAXRED
            KWRK1  = KSOLEQ + MAXRED*MAXRED
            IF (LREDS) THEN
              KREDS = KWRK1 
              KWRK1 = KREDS + MAXRED*MAXRED
            END IF
            LWRK1  = LWORK  - KWRK1
C
            IF (LWRK1.LT. 0 )
     *         CALL QUIT(' TOO LITTLE WORKSPACE IN CC_SOLRSP')

C-------------------------------------------------------------------
C           Frequencies for diagonal scaling 
C-------------------------------------------------------------------
            IF (ISIDE.EQ.1) THEN
              SIGNF = ONE
            ELSE IF (ISIDE.EQ.-1) THEN
              SIGNF = -ONE
            ELSE
              CALL QUIT('Illegal ISIDE parameter in CC_SOLRSP')
            END IF

            IF (LIST(1:2).EQ.'M1' .OR. LIST(1:2).EQ.'N2') THEN
               SIGNE = -ONE
            ELSE
               SIGNE = ONE
            END IF

            DO IFR = 1, NSIMR
               IR  = IRST + IFR - 1
               WORK(KEIVAL+IFR-1) = ZERO
               DO I = 1, NSTAT
                 WORK(KEIVAL+IFR-1) =
     &                 WORK(KEIVAL+IFR-1) + SIGNE * EIGVAL(IR,I)
               END DO
               IF ( LIST(1:1).NE.'C' .AND. LIST(2:2).NE.'C' ) THEN
                 DO I = 1, ORDER
                   WORK(KEIVAL+IFR-1) =
     &                   WORK(KEIVAL+IFR-1) + SIGNF * FREQS(IR,I)
                 END DO
               END IF
            END DO
C
C----------------------
C           Open files.
C----------------------
C
            CALL CC_FILOP(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12,
     *                    FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12,
     *                    FS12AM,LUFS12,FS2AM,LUFS2)
C-------------------------------------------------------------------
C           check for projection: if we solve projected equations
C                                 set flat LPROJECT and store excited
C                                 state index in ISTATPRJ
C-------------------------------------------------------------------
            LPROJECT = .FALSE.
            ISTATPRJ = 0
            IF (LIST(1:2).EQ.'ER' .OR. LIST(1:2).EQ.'EL'
     &                             .OR. LIST(1:3).EQ.'PL1') THEN
              LPROJECT = LPROJ(IRST)
              ISTATPRJ = ISTAT(IRST,1)
              DO IDX = IRST+1,IRST+NSIMR-1
                IF ((LPROJ(IDX).NEQV.LPROJECT) .OR.
     &              (LPROJECT .AND. ISTAT(IDX,1).NE.ISTATPRJ) )
     &           CALL QUIT(
     &           'Error in batching for ERn/ELn/PL1 equations.')
              END DO
            END IF
C
C----------------------------------------------------------------------
C           For Cholesky CC2 we need to delete files for eff. rhs
C           for 1st-order amplitudes (if the files should exist from
C           a previous run).
C----------------------------------------------------------------------
C
            IF ((CHOINT.AND.CC2) .AND. (LIST(1:2).EQ.'R1')) THEN
               ILSTXI = IRST
               DO I = 1,NSIMR
                  CALL CC_RMRSP('eO1',ILSTXI,ISYM)
                  ILSTXI = ILSTXI + 1
               ENDDO
            ENDIF
 
C--------------------------------
C           Create start vectors.
C--------------------------------
C
            CALL CCEQ_STR(FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12,
     *                    FS12AM,LUFS12,FS2AM,LUFS2,LPROJECT,ISTATPRJ,
     *                    TRIPLET,ISIDE,IRST,NSIMR,NUPVEC,
     *                    NREDH,WORK(KEIVAL),WORK(KIPLAC),
     *                    WORK(KWRK1),LWRK1,LIST)
C
            CALL FLSHFO(LUPRI)
C
C---------------------------------------------
C           Solve equations by call to solver.
C---------------------------------------------
C
            IF (LINQCC .AND. CCSDT) ECURR = WORK(KEIVAL)
C
            CALL CCEQ_SOL(LIST,LPROJECT,ISTATPRJ,ECURR,
     *                    FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12,
     *                    FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12,
     *                    LREDS,WORK(KREDS),FS12AM,LUFS12,FS2AM,LUFS2,
     *                    LINQCC,TRIPLET,ISIDE,IRST,NSIMR,NUPVEC,
     *                    NREDH,WORK(KREDH),WORK(KEIVAL),
     *                    WORK(KSOLEQ),WORK(KWRK1),LWRK1,APROXR12)
C
            IF (LINQCC .AND. CCSDT) ECURR = 0.0D0
C
C---------------------------------------------------------
C           Analysis of solution vectors and save on disk.
C---------------------------------------------------------
C
            IF (CHOINT .AND. CC2) THEN

C              Cholesky section.
C              -----------------

               NVARPT = NT1AM(ISYM) + 2*NALLAI(ISYM) + 1

               KWRK2 = KWRK1 + NVARPT
               LWRK2 = LWORK - KWRK2

               IF (LWRK2 .GT. 0) THEN
                  NSIMUL = MIN(NSIMR,LWRK2/NT1AM(ISYM))
                  IF (NSIMUL .LE. 0) THEN
                    CALL QUIT('Insufficient memory in CC_SOLDRV-Chopr2')
                  ENDIF
               ELSE
                  CALL QUIT('Insufficient memory in CC_SOLDRV-Chopr1')
               ENDIF

               NBATCH = (NSIMR - 1)/NSIMUL + 1
               IOFF1  = 1
               ICOUNT = 0
               ILSTNR = IRST

               DO I = 1,NBATCH

                  IOFF2 = MIN(NSIMUL,NSIMR  - (I-1)*NSIMUL)
                  CALL CCCONV(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM,
     *                        TRIPLET,CCR12,NREDH,IOFF1,IOFF2,
     *                        WORK(KSOLEQ),WORK(KWRK2),WORK(KWRK1))

                  DO J = 1,IOFF2

                     ICOUNT = ICOUNT + 1

                     KT1 = KWRK2 + NT1AM(ISYM)*(ICOUNT - 1)

                     IF (IPRINT .GT. 1) THEN
                        WRITE(LUPRI,'(//1X,A)')
     *                  'Analysis of CC Response/Cauchy vector:'
                        WRITE(LUPRI,'(1X,A)')
     *                  '--------------------------------------'
                        CALL CC_PRAM(WORK(KT1),PT1,ISYM,.FALSE.)
                     ENDIF
                     IF ((IPRINT.GT.10).OR.DEBUG) THEN
                        RHO1N = DDOT(NT1AM(ISYM),WORK(KT1),1,
     *                               WORK(KT1),1)
                        WRITE(LUPRI,'(1X,A35,1X,E20.10)')
     *                  'Norm of singles part:   ',RHO1N
                     ENDIF

                     IOPT = 1
                     CALL CC_WRRSP(LIST,ILSTNR,ISYM,IOPT,MODEL,DUMMY,
     &                             WORK(KT1),DUMMY,WORK(KWRK1),NVARPT)
                     IF (LIST(1:2) .EQ. 'R1') THEN
                        CALL CC_RMRSP('eO1',ILSTNR,ISYM)
                     ENDIF
                     ILSTNR = ILSTNR + 1

                  ENDDO

               ENDDO

               GOTO 1234

            END IF
C
C          Conventional section
C
           NVARPT = LETOT + 2*NALLAI(ISYM)
           KWRK2  = KWRK1 + NVARPT
           LWRK2  = LWORK - KWRK2
           NSIMUL = MIN(NSIMR,LWRK2/LETOT )
           NBATCH = (NSIMR -1)/NSIMUL + 1
           IOFF1  = 1
           ICOUNT = 0
           ILSTNR = IRST
C
           DO I = 1,NBATCH
              IOFF2 = MIN(NSIMUL,NSIMR  - (I-1)*NSIMUL)
              CALL CCCONV(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM,
     *                    TRIPLET,CCR12,NREDH,IOFF1,IOFF2,WORK(KSOLEQ),
     *                    WORK(KWRK2),WORK(KWRK1))
              IF (DEBUG) THEN
                 WRITE(LUPRI,*) 'After CCCONV '
                 WRITE(LUPRI,*) 'IOFF1,IOFF2,NREDH',IOFF1,IOFF2,NREDH
                 WRITE(LUPRI,*) 'KWRK1,KWRK2',KWRK1,KWRK2
                 WRITE(LUPRI,*) 'LWRK1,LWRK2',LWRK1,LWRK2
                 X1  = DDOT(LETOT,WORK(KWRK1),1,WORK(KWRK1),1)
                 WRITE(LUPRI,*) 'norm KWRK1:',X1
                 X2  = DDOT(LETOT,WORK(KWRK2),1,WORK(KWRK2),1)
                 WRITE(LUPRI,*) 'norm KWRK2:',X2
              ENDIF
C
              IF ((IPRINT .GT. 30) .OR. DEBUG) THEN
                 CALL AROUND('CC_SOLRSP: response vectors in mo basis')
                 CALL OUTPUT(WORK(KWRK2),1,NCCVAR,1,IOFF2,
     *                       NCCVAR,IOFF2,1,LUPRI)
              ENDIF
C
              DO J = 1,IOFF2
                 ICOUNT = ICOUNT + 1
C
                 KT1 = KWRK2 + (ICOUNT-1)*LETOT
                 KT2 = KT1 + NT1AM(ISYM)
                 IF (CCR12) KT12 = KT2 + NT2AM(ISYM)
                 IF (IPRINT .GT. 1) THEN
                    WRITE(LUPRI,'(//1X,A)')
     *         'Analysis of CC Response/Cauchy vector:'
                    WRITE(LUPRI,'(1X,A)')
     *         '--------------------------------------'
C
                    CALL CC_PRAM(WORK(KWRK2 + (ICOUNT-1)*LETOT),
     *                           PT1,ISYM,.FALSE.)
                    IF ((IPRINT.GT.2).OR.DEBUG) THEN
                      RHO1N = DSQRT(DDOT(NT1AM(ISYM),WORK(KT1),1,
     *                                               WORK(KT1),1))
                      WRITE(LUPRI,'(1X,A35,1X,E20.10)')
     *               'Norm of singles part:   ',RHO1N
                      IF ( .NOT. CCS) THEN
                       RHO2N = DSQRT(DDOT(NT2AM(ISYM),WORK(KT2),1,
     *                                                WORK(KT2),1))
                       WRITE(LUPRI,'(1X,A35,1X,E20.10)')
     *                 'Norm of doubles part:   ',RHO2N
                      ENDIF
                      IF (CCR12) THEN
                       RHO12N = DSQRT(DDOT(NTR12AM(ISYM),WORK(KT12),1,
     &                                                   WORK(KT12),1))
                       WRITE(LUPRI,'(1X,A35,1X,E20.10)')
     *                 'Norm of r12 doubles part:   ',RHO12N
                      ENDIF
                    ENDIF
                 ENDIF
C
C----------------------------------------------
C                Save response vectors on file.
C----------------------------------------------
C
                 KT1    = KWRK2 + (ICOUNT-1)*LETOT
                 KT2    = KT1 + NT1AM(ISYMTR)

                 IF (CCSTST) THEN
                   CALL DZERO(WORK(KT2),NT2AM(ISYMTR))
                 END IF

                 IOPT   = 3
                 CALL CC_WRRSP(LIST,ILSTNR,ISYM,IOPT,MODEL,DUMMY,
     *                  WORK(KT1),WORK(KT2),WORK(KWRK1),NVARPT)
                 IF (CCR12) THEN 
                   IOPT = 32
                   KT12 = KT2 + NT2AM(ISYMTR)
                   CALL CC_WRRSP(LIST,ILSTNR,ISYM,IOPT,MODEL,DUMMY,
     *                    DUMMY,WORK(KT12),WORK(KWRK1),NVARPT)
                 END IF
                 ILSTNR = ILSTNR + 1
C
              END DO
              IOFF1 = IOFF1 + NSIMUL
           END DO
C
 1234      CONTINUE       ! Cholesky
C
C---------------------------------
C          Close and delete files.
C---------------------------------
C
           CALL CC_FILCL(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12,
     *                   FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12,
     *                   FS12AM,LUFS12,FS2AM,LUFS2)
C
         END DO


      ELSE
        WRITE(LUPRI,*) 'There are no amplitudes of symmetry ',ISYM,'!'
      ENDIF
C
      END DO

C---------------------------------------------------------------------
C     For CC3 restore NSIMLE
C---------------------------------------------------------------------
      IF (CCSDT) NSIMLE = NSIMLEBAK

C-------------------------------------------------------------
C     finished
C-------------------------------------------------------------
      CALL QEXIT('CC_SOLDRV')
      RETURN
      END
C--------------------------------------------------------------------
C  /* Deck cc_getvec */
      SUBROUTINE CC_GETVEC(LUC1,FC1AM,LUC2,FC2AM,LUC12,FC12AM,
     *                     TRIPLET,IVEC,VEC)
C
C     Ove Christiansen 5-11-1996:
C
C     Get vector nr. IVEC from files 
C            FC1AM  (LUC1)  singles part
C            FC2AM  (LUC2)  doubles part
C            FC12AM (LUC12) R12 doubles part
C
      IMPLICIT NONE
C#include "implicit.h"
#include "priunit.h"
C#include "iratdef.h"
#include "cclr.h"
#include "ccsdinp.h"
#include "ccsdsym.h"
Cholesky
#include "maxorb.h"
#include "ccdeco.h"
Cholesky

      DOUBLE PRECISION VEC(*), DDOT

      LOGICAL TRIPLET
      CHARACTER FC1AM*(*), FC2AM*(*), FC12AM*(*)
      INTEGER LUC1, LUC2, LUC12, IVEC, NT2VAR, NT12VAR, IADR
Cholesky
      LOGICAL CCSEFF
C
      CCSEFF = CCS .OR. (CHOINT .AND. CC2)
Cholesky
C
      CALL CC_RVEC(LUC1,FC1AM,NT1AM(ISYMTR),NT1AM(ISYMTR),IVEC,VEC(1))
Chol  IF (.NOT. CCS) THEN
      IF (.NOT. CCSEFF) THEN
         IF (TRIPLET) THEN
            NT2VAR = NT2AMT(ISYMTR)
         ELSE
            NT2VAR = NT2AM(ISYMTR)
         END IF                          
         IADR = 1 + NT1AM(ISYMTR)
         CALL CC_RVEC(LUC2,FC2AM,NT2VAR,NT2VAR,IVEC,VEC(IADR))
         IF (CCR12) THEN
            NT12VAR = NTR12AM(ISYMTR)
            IADR    = 1 + NT1AM(ISYMTR) + NT2VAR
            CALL CC_RVEC(LUC12,FC12AM,NT12VAR,NT12VAR,IVEC,VEC(IADR))
         END IF
         IF (CCSTST) 
     *     CALL DZERO(VEC(1+NT1AM(ISYMTR)),NCCVAR-NT1AM(ISYMTR))
      ENDIF

      RETURN
      END
C  /* Deck cc_putvec */
      SUBROUTINE CC_PUTVEC(LUC1,FC1AM,LUC2,FC2AM,LUC12,FC12AM,
     *                     TRIPLET,IVEC,VEC)
C
C     Ove Christiansen 22-1-1996:
C
C            put vector nr. IVEC to 
C                singles file (name: FC1AM,  nr: LUC1)
C                doubles file (name: FC2AM,  nr: LUC2)
C            R12 doubles file (name: FC12AM, nr: LUC12)
C
      IMPLICIT NONE
C#include "implicit.h"
#include "priunit.h"
C#include "iratdef.h"
#include "cclr.h"
#include "ccsdinp.h"
#include "ccsdsym.h"
Cholesky
#include "maxorb.h"
#include "ccdeco.h"
Cholesky

      DOUBLE PRECISION VEC(*), DDOT

      LOGICAL TRIPLET
      CHARACTER FC1AM*(*),FC2AM*(*),FC12AM*(*)
      INTEGER LUC1, LUC2, LUC12, IVEC, NT2VAR, NT12VAR, IADR
C
Cholesky
      LOGICAL CCSEFF
C
      CCSEFF = CCS .OR. (CHOINT .AND. CC2)
Cholesky
C
      CALL CC_WVEC(LUC1,FC1AM,NT1AM(ISYMTR),NT1AM(ISYMTR),
     *             IVEC,VEC(1))
Chol  IF ( .NOT. CCS) THEN
      IF ( .NOT. CCSEFF) THEN
         IF (CCSTST) 
     *     CALL DZERO(VEC(1+NT1AM(ISYMTR)),NCCVAR-NT1AM(ISYMTR))
         IF (TRIPLET) THEN
            NT2VAR = NT2AMT(ISYMTR)
         ELSE
            NT2VAR = NT2AM(ISYMTR)
         END IF                              
         IADR = 1 + NT1AM(ISYMTR)
         CALL CC_WVEC(LUC2,FC2AM,NT2VAR,NT2VAR,IVEC,VEC(IADR))
         IF (CCR12) THEN
            NT12VAR = NTR12AM(ISYMTR)
            IADR    = 1 + NT1AM(ISYMTR) + NT2VAR
            CALL CC_WVEC(LUC12,FC12AM,NT12VAR,NT12VAR,IVEC,VEC(IADR))
         END IF
      ENDIF
C
      RETURN
      END
C  /* Deck cc_filop  */
      SUBROUTINE CC_FILOP(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12,
     *                    FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12,
     *                    FS12AM,LUFS12,FS2AM,LUFS2)
C
C     Ove Christiansen 22-1-1996:
C                Open files for linear solver.
C                Has to be 8 character long.
C
      IMPLICIT NONE
C#include "implicit.h"
C#include "iratdef.h"
#include "ccsdinp.h"
Cholesky
#include "maxorb.h"
#include "ccdeco.h"
#include "priunit.h"
#include "ccsdsym.h"
#include "ccorb.h"
      LOGICAL CCSEFF, OLDDX
      INTEGER MAXT1, ISYM, LRLEN
C
      INTEGER MAX
Cholesky
      INTEGER LUFR1,LUFR2,LUFR12,LUFC1,LUFC2,LUFC12,LUFS12,LUFS2
      CHARACTER*8 FRHO1,FRHO2,FC1AM,FC2AM,FRHO12,FC12AM,FS12AM,FS2AM
C
      IF (CHOINT) THEN

C        Cholesky (fortran) i/o.
C        Use max. T1 dimension as record length.
C        ---------------------------------------

         MAXT1 = NT1AM(1)
         DO ISYM = 2,NSYM
            MAXT1 = MAX(MAXT1,NT1AM(ISYM))
         ENDDO
         LRLEN = 2*MAXT1
         LRLEN = MAX(LRLEN,1)

         CALL GPOPEN(LUFC1,FC1AM,'UNKNOWN','DIRECT','UNFORMATTED',
     &               LRLEN,OLDDX)
         CALL GPOPEN(LUFR1,FRHO1,'UNKNOWN','DIRECT','UNFORMATTED',
     &               LRLEN,OLDDX)

      ELSE

C        Standard (cray) i/o.
C        --------------------
C
         CCSEFF = CCS .OR. (CHOINT.AND.CC2) ! kept if we should switch later

         CALL WOPEN2(LUFC1,FC1AM,64,0)
         CALL WOPEN2(LUFR1,FRHO1,64,0)
C
Chol     IF ( .NOT. CCS ) THEN
         IF ( .NOT. CCSEFF ) THEN
            CALL WOPEN2(LUFC2,FC2AM,64,0)
            CALL WOPEN2(LUFR2,FRHO2,64,0)
         ENDIF
C
         IF ( CCR12 ) THEN
            CALL WOPEN2(LUFR12,FRHO12,64,0)
            CALL WOPEN2(LUFC12,FC12AM,64,0)
            CALL WOPEN2(LUFS12,FS12AM,64,0)
            CALL WOPEN2(LUFS2,FS2AM,64,0)
         END IF
C
      END IF
C
      RETURN
      END
C  /* Deck cc_filcl  */
      SUBROUTINE CC_FILCL(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12,
     *                    FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12,
     *                    FS12AM,LUFS12,FS2AM,LUFS2)
C
C     Ove Christiansen 22-1-1996:
C
C     tbp: use fortran i/o [notably, gpopen/gpclose] for
C          Cholesky....
C
      IMPLICIT NONE
C#include "implicit.h"
C#include "iratdef.h"
#include "ccsdinp.h"
Cholesky
#include "maxorb.h"
#include "ccdeco.h"
      LOGICAL CCSEFF
Cholesky
C
      INTEGER LUFR1,LUFR2,LUFR12,LUFC1,LUFC2,LUFC12,LUFS12,LUFS2
      CHARACTER*8 FRHO1,FRHO2,FC1AM,FC2AM,FRHO12,FC12AM,FS12AM,FS2AM
C
      IF (CHOINT) THEN

C        Close files (fortran).
C        ----------------------

         CALL GPCLOSE(LUFC1,'DELETE')
         CALL GPCLOSE(LUFR1,'DELETE')

      ELSE
C
C-----------------------------
C        Close files. (crayio)
C-----------------------------
C
         CCSEFF = CCS .OR. (CHOINT.AND.CC2) ! kept if we should switch later
C
         CALL WCLOSE2(LUFC1,FC1AM,'DELETE')
         CALL WCLOSE2(LUFR1,FRHO1,'DELETE')
C
Chol     IF ( .NOT. CCS) THEN
         IF ( .NOT. CCSEFF) THEN
            IF (LUFC2 .GT. 0) CALL WCLOSE2(LUFC2,FC2AM,'DELETE')
            IF (LUFR2 .GT. 0) CALL WCLOSE2(LUFR2,FRHO2,'DELETE')
         ENDIF
C
         IF ( CCR12 ) THEN
            CALL WCLOSE2(LUFR12,FRHO12,'DELETE')
            CALL WCLOSE2(LUFC12,FC12AM,'DELETE')
            CALL WCLOSE2(LUFS12,FS12AM,'DELETE')
            CALL WCLOSE2(LUFS2,FS2AM,'DELETE')
         END IF
C
      END IF
C
      RETURN
      END
C--------------------------------------------------------------------
C
C  /* Deck cc_chonsim */
      SUBROUTINE CC_CHONSIM(LIST,ISYM,NEQSYM,NSIMLE,NSIM)
C
C     Thomas Bondo Pedersen, April 2003.
C
C     Purpose: Set the number of equations to be solved simultaneously
C              for Cholesky CC2.
C
#include "implicit.h"
      CHARACTER*(*) LIST
#include "cclrcho.h"
#include "priunit.h"

      CHARACTER*10 SECNAM
      PARAMETER (SECNAM = 'CC_CHONSIM')

      IF (LIST(1:2) .EQ. 'L0') THEN

         NSIM = 1

      ELSE IF (LIST(1:2) .EQ. 'R1') THEN

         IF (NSIMLE .EQ. 0) THEN
            NSIM = NCHOPLR(ISYM)
         ELSE
            NSIM = NSIMLE
            IF ((NSIM.NE.1) .AND. (NSIM.NE.NCHOPLR(ISYM))) THEN
               WRITE(LUPRI,'(/,1X,A,A,/)')
     &         '*** NOTICE: .NSIMLE input ignored - ',
     &         'Cholesky R1 equations solved 1 at a time!'
               NSIM = 1
            ENDIF
         ENDIF

      ELSE

         WRITE(LUPRI,'(//,A,A,A,/)')
     &   LIST(1:3),' not recognized in ',SECNAM
         CALL QUIT('Unknown list in '//SECNAM)

      ENDIF

      IF ((NSIM.LE.0) .OR. (NSIM.GT.NEQSYM)) THEN
         WRITE(LUPRI,*)
         WRITE(LUPRI,*) 'Error in ',SECNAM,':'
         WRITE(LUPRI,*) 'NSIMLE  = ',NSIMLE
         WRITE(LUPRI,*) 'NEQSYM  = ',NEQSYM
         WRITE(LUPRI,*) 'NCHOPLR = ',NCHOPLR(ISYM),' (ISYM = ',ISYM,')'
         WRITE(LUPRI,*) 'NSIM    = ',NSIM
         WRITE(LUPRI,*) '(This is a Cholesky-specific problem)'
         CALL QUIT('NSIM error in '//SECNAM)
      ENDIF

      RETURN
      END
